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NOMENCLATURE 


A,  B,  C 

a 

h 

c 

d 

f 

A' 


Rs,  gi 
h 
hi 
hi 

K 

Ko 

k 

P 

T 

T, 

Ts 

t 

u,  V 
U*,  V* 

u',  v' 


constants 

combined  convection-diffusion  coefficient 
source  term  in  the  discretization  equation 
specific  heat 
coefficient 
mass  fraction 

acceleration  due  to  gravity 

volume  fraction  of  the  solid  phase  and  liquid  phase,  respectively 

enthalpy 

solidus  enthalpy 

iiquidus  enthalpy 

heat  of  fusion 

permeability 

constant  for  permeability 

thermal  conductivity 

pressure 

temperature 

Iiquidus  temperature 

solidus  temperature 

time 

velocity  components 
guessed  velocity  components 
velocity  corrections 


vii 


u,  V  pseudo- velocity  components 

y  velocity  vector 

ji,  y  Cartesian  coordinates 

Greek  Symbols 

P  thermal  expansion  coefficient 

fi  dynamic  viscosity 

p  density 

p  partial  density 

Subscripts 

E,  IV,  N,  S  cast,  west,  north,  south  neighbors  of  the  control  volume  center 
e,  w,  ft,  s  cast,  west,  north,  and  south  faces  of  a  control  volume 

/  liquid  phase 

P  pertaining  to  the  grid  point,  P 

p  pressure 

5  solid  phase 

u,  V  pertaining  to  the  velocity  components,  u  and  v 

X,  y  pertaining  to  the  Cartesian  coordinates 
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I.  INTRODUCTION 


Processes  related  to  phase  change  encompass  a  wide  range  of  engineering 
and  scientific  disciplines  and  occur  in  many  applications  such  as  welding,  cast¬ 
ing  and  energy  storage.  Owing  to  the  absorption  or  release  of  latent  thermal 
energy,  phase  change  problems  arc  nonlinear,  and  exact  solutions  arc  limited  to 
a  small  class  of  problems  involving  pure  substances  in  one-dimensional  infinite 
or  semi-infinite  domains  [1,2].  Unfortunately,  most  pr'.ctical  phase  change 
problems  are  multi-dimensional;  the  thermophysical  properties  such  as  thermal 
conductivity,  density,  and  specific  heat  for  the  phase  change  material  (PCM)  are 
changing  with  the  changing  of  phase,  and  free  convection  occurs  in  melting  liq¬ 
uid  PCM.  These  have  focused  attention  on  development  of  suitable  numerical 
procedures. 

Generally,  the  numerical  techniques  can  be  divided  into  two  groups.  The 
first  group  uses  the  front  tracking  method,  which  utilizes  two  independent  con¬ 
servation  equations  for  each  phase  and  couples  them  with  appropriate  boundary 
conditions  at  the  phase  interface.  Front  tracking  methods  require  the  existence 
of  discrete  interfaces  between  phases  in  the  domain  and  are  generally  limited  to 
pure  substances.  The  primary  difficulty  associated  with  this  method  centers  on 
tracking  the  phase  interface,  which  is  generally  an  unknown  function  of  space 
and  time.  The  need  for  moving  numerical  grids  and/or  coordinate  mapping 
procedures  also  complicates  the  application  of  this  technique. 

The  second  group  uses  the  enthalpy  method  which  is  sometimes  called  the 
fixed  grid  method.  The  enthalpy  method  does  not  track  the  phase  interface, 
instead,  simply  calculates  the  enthalpy  of  the  PCM  at  each  numerical  grid. 
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Only  one  set  of  conservation  equations  arc  needed  for  both  solid  and  liquid 
PCM  domains.  This  eliminates  the  complications  of  tracking  phase  interface 
and  moving  numerical  grids.  Moreover,  the  enthalpy  method  allows  a  mushy 
region  (solid-liquid  phases  coexist  region).  This  enables  the  enthalpy  method  to 
be  implemented  for  multiconstitucnt  systems,  which  do  not  exhibit  a  sharp  in¬ 
terface  between  the  solid  and  liquid  phases.  The  enthalpy  method  is  generally 
easier  to  be  implemented  for  solidification  of  alloys,  melting  of  impure  sub¬ 
stances. 

Based  on  enthalpy  formulation,  several  models  were  developed  to  solve 
phase  change  problems.  The  first  one  is  the  conduction  model  [3,4].  This  is  the 
easiest  model,  where  free  convection  effects  of  the  liquid  PCM  are  not  consid¬ 
ered.  With  the  exception  of  microgravity  applications,  it  is  frequently  this  frcc- 
convcctivc  motion  which  is  the  dominant  mode  of  heat  transfer.  To  include  free 
convection  influences,  Schneider  [5]  proposed  a  numerical  model  to  solve  phase 
change  problems  for  pure  substances.  By  assuming  the  behavior  of  fluid  flows 
in  the  mushy  region  to  be  similar  to  that  of  fluid  flows  in  porous  media,  Voller 
ct  al.  [6,7]  proposed  the  Enthalpy-Porosity  model.  This  model  was  developed 
by  utilizing  Darcy's  law  to  modify  the  pressure  gradients  in  the  momentum 
equations.  A  continuum  model  for  analyzing  the  solid-liquid  phase  change  be¬ 
havior  in  a  binary  system  was  developed  by  Bennon  and  Incropera  [8  —  10]. 
Semi-empirical  laws,  as  well  as  microscopic  descriptions  of  the  transport  be¬ 
havior,  have  been  integrated  with  the  principles  of  classical  mixture  theory  in 
obtaining  this  model. 
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In  the  present  research,  a  continuum  model  is  utilized  and  a  finite-differ¬ 
ence  numerical  method  is  used  to  solve  the  conservation  equations. 

A.  CHANGING  OF  PHYSICAL  PROPERTIES 

Physical  properties  of  the  PCM  such  as  density,  specific  heat,  thermal 
conductivity,  etc.  change  while  undergoing  solid-liquid  phase  change.  In  this 
research,  a  numerical  method  will  be  developed  to  accommodate  the  changing 
of  thc.se  physical  properties. 

B.  SLUMPING  FLOATING  PHENOMENON 

The  variation  of  solid  and  liquid  densities  can  cause  slumping/floating 
phenomenon  during  the  melting  processes.  Floating  of  ice  in  the  water  is  an 
example.  Slumping  or  floating  of  solid  PCMs  in  enclosures  can  significantly  en¬ 
hance  the  melting  rate  [11,12].  Simple  analyses  were  developed  by  presuming 
the  shapes  of  solid  PCM  [13,14],  or  by  considering  that  conduction  is  the  only 
energy  transport  mechanism  [15].  However,  no  numerical  studies  for  melting 
processes  including  both  free  convection  and  slumping  have  been  reported.  The 
difficulties  are  due  to  the  fact  that  moving  of  the  solid  domain  is  unknown  and 
time  dependent.  The  utilizing  of  a  continuum  model  enables  the  prediction  of 
the  velocity  for  each  phase,  and  provides  the  possibility  of  solving  the  melting 
heat  transfer  problem  having  a  slumping/floating  phenomena. 


C.  DENSITY  VARIATION  INDUCED  FLOW 


During  a  melting  process,  the  variation  of  phase  densities  can  cause  motion 
of  the  liquid  PCM.  If  the  density  of  the  liquid  PCM  is  smaller  than  that  of  the 
solid  PCM,  an  expansion-induced  flow  would  occur.  Variation  of  phase  densi¬ 
ties  is  common,  and  for  many  materials,  the  variations  are  large.  For  example, 
the  liquid  density  of  PI  16  wax  is  760  kglm^,  and  the  solid  density  is  818 
However,  the  importance  of  density  variation  induced  flow  during  melting  has 
not  yet  been  reported.  In  this  research,  the  effect  of  density  variation  induced 
flow  will  be  studied. 

Experimental  results  will  be  used  to  verify  the  reliability  of  the  numerical 
solutions.  PI  16  wax  with  a  melting  temperature  of  46.7  C  (1 16  F)  and  a  mushy 
temperature  range  of  10  C  [16]  has  been  selected  for  the  experiments.  Physical 
properties  of  this  wax  will  be  used  in  the  numerical  calculations. 
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II.  MATHEMATICAL  FORMULATION 

A.  CONSERVATION  OF  MASS 

The  principle  of  conservation  of  mass  is  that  the  time-rate  of  change  of 
mass  inside  a  control  volume  equals  the  net  integral  of  the  mass  flux  over  the 
control  volume.  The  two-dimensional,  solid-liquid  system,  according  to  the 
principle  of  conservation  of  mass,  yields 

+  P/ “/]  ^feP,v,  +g/P/V/]  0  tn 

dt  dx  dy  ^  ^ 

This  can  be  seen  from  Fig.  1.  Let  the  mixture  density,  p,  be  defined  by 

P  =  SsPs  +^/P/ 

=  Ps  +  Pi  (  2  ) 


and  the  mass  averaged  velocity,  F,  be  defined  by 


V  = 


-  fs  K  +fi  ^I 


(3) 


then  substituting  Eqs.(2)  and  (3)  into  Eq.(l),  it  follows  that 


dp  d{pu) 

dt  dx 


d{pv) 

dy 


=  0 


(4) 
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This  is  the  equation  describing  the  conservation  of  mass  for  the  two-dimen 
sional,  solid-liquid  system. 


Vs  Vl 


Figure  1.  A  control  volume  in  mushy  region. 


B.  CONSERVATION  OF  MOMENTUM 

The  equations  describing  the  conservation  of  momentum  for  the  two-di' 
mcnsional  system  are 


dt 


(5 


where  K  is  the  permeability.  Details  of  derivation  of  these  equations  are  shown 
in  the  Appendix.  In  the  present  analysis,  the  permeability  is  assumed  to  vary 
with  liquid  volume  fraction  according  to  the  Kozeny-Carman  equation  [17] 


K  =  —  ~  ]  (7) 

(1 

where  Kq  is  a  constant  which  depends  on  the  specific  multiphase  region  mor¬ 
phology. 

C.  CONSERVATION  OF  ENERGY 

The  equation  describing  the  conservation  of  energy  is 

+  V.  =  V.  (-f  V/t) -t- V.  [-^  - /!)] 

dt  c,  c, 

where  the  mixture  enthalpy  and  thermal  conductivity  are 


h=fA-^fihi 

(9) 

^  +  Siki 

(10) 

Phase  enthalpies  are 


^5  = 


cj 


(11) 


hi  =  CiT  -I-  [(c,  -  c/)r,  +  hj\ 


(12) 
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where  it  is  presumed  that  hs\  r=o^O  and  that  (/?/ - |7=7- = ///.  Details  of 
derivation  of  the  energy  eciualion  arc  shown  in  the  Appenilix. 

D.  THERMODYNAMIC  RELATIONS 


By  assuming  that  the  phase  fractions  arc  linear  functions  of  temperature  in 
the  mushy  range  and  in  the  saturated  condition,  we  can  obtain 

r  -  r 

/5=>-//  (14) 

From  Eq.(2),  w'c  can  find  the  phase  volume  fractions, 

P- Pi 

Ss  p^-pi  (*5) 

(16) 

Let  h[  and  /i2  represent  the  enthalpies  of  the  PCM  at  the  temperatures, 
Tv  and  T/,  respectively.  By  definitions, 

(17) 


and 


^2  =  +  Ci{T,  -  T)  +  hf 


(18) 


8 


For  h<h\,  wc  can  determine  the  temperature  of  the  PCM  by  using  the  follow¬ 
ing: 


T  = 


(19) 


For  /?|  ^  /?  <  hi,  from  Eqs.(9),  (11-14),  we  can  determine  the  temperature  of  the 
PCM  by  using  the  following: 


AT'^  BT C  =  0 


where 


A 

B 

C 


Cs  -  Cl 
T,-Ts 

-lA{Ti+Ts)  + 

A 


Ti- Ts 


+  C/] 


Ti-Ts 


Ts  +  h 


(20) 


For  h>h2,  from  Eq.(12),  we  can  determine  the  temperature  of  the  PCM  by 
using  the  following: 


(21, 
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III.  NUMERICAL  METHOD 


A.  STAGGERED  GRID 

Figure  2  shows  staggered  locations  for  u,  v,  and  p.  A  staggered  control 
volume  for  the  x-momentum  equation  is  shown  in  Fig  3.  Figure  4  shows  the 
control  volume  for  the  y-momentum  equation. 


Figure  2.  Staggered  locations  for  computational  variables. 


Figure  3.  Control  volume  for  x-momentum  equation 


Figure  4.  Control  volume  for  y-momentum  equation 


B.  NUMERICAL  FORMULATION 


The  discretized  equations  for  the  continuity,  x-momentum,  and  y-momcn- 
tum  equations  are 


[pp  - 


+  ^yVPe^ip  -  +  ^x(p„vp  -  p^v^)  =  0 


-  ^uE^E  +  ^uW^W  T  ^u.v“,V  +  ^uS^S  +  +  A>'(pp  -  Pp) 

=  '^^unb^nb  +  ^x  +  ^yiPp  -  Pe) 

u^.pyp  =  «v/,v/,-  +  +  ^v.vLv  +  ^vs'^s  +  +  ^iPp  -  P{^) 

=  Y^<^vnb^nb  +  +  ^'^iPp  ”  Ps) 

Let  the  pseudo-\clocitics  u  and  v  be  defined  by  the  following; 


Yj^unb^, 


nb  ^  ^x 


'y  '/^vnb^nb  T  ^y 


and  substituting  into  Eqs.(23)  and  (24),  it  follows  that 


Up  =  U  +  d^pipp  -  pi,) 


'>  -  ^'n  +  ^piPp  ~  P\) 
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where 


d^,,  -  ^y|a^p 


d,p  =  Axja^p. 


i 


Substituting  Eqs.(27)  and  (28)  into  Eq.(22),  it  follows  that 


^pPPp  -  ^^pnbPnb  +  ^P 


(29) 


where 


{pp-  Pp)  Ax  Ay 
At 


+  Ay(P  A  -  pjiv)  +  ^{Pn^p  -  Psis)'l 


(30) 


Let  u*,  V*  and  p*  be  the  guessed  values  and 


^uP^*P  =  Yj^unb^*nb  +  +  ^y(P*P  -  P*e)  •  ) 

^.P^*P  =  +  by  +  Ax{p*p  -  p*f^)  (32) 

It  is  proposed  that 


-  ^piP'p  ~  P'e) 

(33) 

=  ^pip'p  ~  P'n) 

(34) 
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where  u',  v',  and  p'  are  corrections  for  velocity  and  pressure,  respectively,  and 
the  correct  velocities  are  obtained  from 


u  —  u*  +  u'  (35) 

V  =  V*  +  v'  (36) 

Substituting  Eqs.(33)  and  (34)  into  Eqs.(35)  and  (36),  respectively,  it  follows 
that 

Up  —  u*  p  +  d^pip' P  —  p' (37) 
V  p  =  V*  p  +  d^pip' p  —  p' i^)  (3S) 

Substituting  Eqs.(37)  and  (38)  into  Eq.(22),  it  follows  that 

^pPP'p  =  YfpnhP'nb  +  ^P  (39) 


where 


[Pp- 

At 


+  Ay{PeU*p  -  p^u*,y)  +  Ax{p„v*p  -  p^v*^)] 


(40) 


C.  NUMERICAL  PROCEDURES 


The  SIMPLER  (Scmi-lmplicit  Method  for  Pressure- Linked  Equations  - 
Revised)  [18]  scheme  is  proposed  to  solve  the  conservation  equations.  The  nu¬ 
merical  procedures  are 
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1.  Given  initial  and  boundary  conditions. 

2.  Advance  a  time  step. 

3.  Given  initial  calculation  values  for  velocities. 

4.  Calculate  pseudo-velocities  u  and  v. 

5.  Calculate  pressure  field. 

6.  Use  the  newest  pressure  as  guessed  pressure  to  find  u*  and  v*. 

7.  Calculate  pressure  corrections  p'. 

8.  Obtain  velocities  u  and  v. 

9.  Find  solutions  for  energy  equation. 

10.  Return  to  step  4  and  repeat  until  convergence. 

1 1 .  Return  to  step  2  for  next  time  step  .solutions  or  stop. 
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IV.  RESULTS 


Figure  5  is  an  example  of  the  numerical  results  showing  the  velocity  fields 
and  the  liquidus  lines.  PI  16  wax  is  the  PCM  constrained  in  a  2.5-cm  x  2.5-cm 
square  enclosure.  The  initial  temperature  's  35  C,  and  faces  are  subjected  to 
constant  temperature  (50  C)  heating. 


1200  seconds  1570  seconds 


Figure  5.  Velocity  vectors  and  liquidus  lines  during  melting. 
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V.  SUMMARY 


I  he  icseareh  is  summarized  below: 

This  research  is  aimed  at  numerical  examinations  of  melting  heat  transfer 
in  enclosures. 

A  continuum  model  is  used  to  describe  the  conservation  equations  for 
the  PCM  during  melting. 

SIMPLER,  a  control-volume-based  numerical  scheme,  is  utilized  to  solve 
the  coupled  conservation  equations. 

Influence  of  the  density  variation  between  solid  and  liquid  phases  on  the 
melting  process  will  be  studied. 

Influence  of  the  slumping/floating  phenomena  on  the  melting  process 
will  be  investigated. 

Solid-liquid  interfacial  motion  experiments  will  be  performed. 
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APPENDIX 


DERIVATION  OF  CONSERVATION  EQUATIONS 


A.  CONSERVATION  OF  MASS 


The  principle  of  conservation  of  mass  states  that  the  time  rate  of  change 
of  mass  for  a  control  volume  equals  the  net  mass  flux  in  the  control  volume. 
For  a  two-dimensional,  solid-liquid  system,  the  conservation  of  mass  equation 
is 

'TA'./'.  +  +/>'// VW/]  ,  n 

- ^ -  +  - r -  -h  - -  =  U  I  I  j 

(Ir  dx  dy  ' 

The  mass  density  of  the  mixture,  p,  is  defined  by 

P  =  SsPs  +S1P1 


-  Ps  Pi 

and  the  mass  averaged  velocity,  F,  is  defined  by 


(2) 


=  fs  +fl  ^1 


Substituting  Eqs.(2)  and  (3)  into  Eq.(l),  we  obtain 


(3) 
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dt  ^  dx 


(4) 


Thus,  this  equation  describes  the  conservation  of  mass  for  the  two-dimensional, 
solid-liquid  system.  In  the  liquid  region,  the  form  of  this  equation  is  the  same 
as  the  traditional  continuity  equation.  If  phase  densities  are  constants,  the  term. 
dpjet,  would  be  nonzero  only  in  the  mushy  zone.  However,  in  the  mushy  zone, 
if  the  solid  density  and  the  liquid  density  are  different  and  if  the  solid  and  liquid 
fractions  arc  time  dependent,  then  the  mass  density  of  the  mixture  is  time  de¬ 
pendent.  The  flow  from  the  density  variation  in  the  solid-liquid  phase  change 
is  due  mainly  to  the  variation  of  solid  and  liquid  densities.  For  example,  during 
the  solidification  of  many  metals,  shrinkage  induced  flow  can  result. 


B.  CONSERVATION  OF  MOMENTUM 


The  momentum  equations  arc  derived  from  Newton's  Second  Law,  which 
states  that  the  product  of  mass  and  acceleration  is  equal  to  the  sum  of  the  ex¬ 
ternal  forces  acting  on  the  body.  Thus,  the  x-momentum  equation  of  the  two- 
dimensional,  solid-liquid  system  is 


{pi  Hi  +  Ps  =  V  .  (g/  au  +  gs  Osx)  +  {pi  Bu  +  Ps  Bsx) 

+  iSi  ^Ix  +  Ss  <^sx) 

I 


(5) 


22 


The  flux  vector,  ax,  represents  the  component  of  the  general  material  stress 
tensor  which  influences  the  x-direction  momentum,  while  Bx  represents  the  x- 
component  body  force  and  Gx  is  the  momentum  production  owing  to  phase  in¬ 
teractions. 

The  left-hand  side  of  Eq.(5)  represents  the  product  of  mass  and  acceler¬ 
ation.  It  can  be  decomposed  as  follows: 

n  d(pu) 

(P/ +  Ps  ^s)  =  +  V  ■  (P/  «/  +  Ps  ^s)  (  ^  ) 

Substituting  Eq.(6)  into  Eq.(5),  it  follows  that 

+  V  .  (p/  Vi  u,  +  i;  =  V  .  (g,  a^^)  +  p  (  7  ) 

where 

B.  =  -^B,,  +  .^B„  =  /,B„  +/,B„  (8) 

and 

=  Si  Gix  +  Ss^sx  (  9  ) 

The  advective  momentum  flux  can  be  decomposed  as  follows: 

(p/  E/  ui  +  Ps  Vs  Us)  =  p  V  u  ->i  ps{Vs  -  V)  (us  -  u) 

+  P/(f''/  -  (10) 

Substituting  Eq.(lO)  into  Eq.(7),  it  follows  that 
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+  V  -  (/J  I  'u)  --  V  -  (^'/  O/x  +  Jk'v  <7.vjt) 


-  V  •  n{u^  ~  u)  +  p,{V,~  V){ui-u)-]  +  +F,  (  II  ) 

The  x-component  of  the  material  stress  tensor  can  be  separated  into  isotropic 
and  deviatoric  components, 


=  -pi  +  ('2) 

It  is  important  to  recognize  that  includes  only  stresses  resulting  from  inter¬ 
action  of  a  single  phase  with  itself.  The  effect  of  interactions  between  phases  is 
accommodated  by  the  quality,  Fx-  Specification  of  requires  a  priori  assess¬ 
ment  of  the  continuity  of  each  mixture  phase.  A  phase  is  considered  to  be  con¬ 
tinuous  if  any  two  points  within  the  phase  can  be  joined  by  a  continuous  curve 
which  lies  solely  within  the  phase.  If  each  pha.se  mixture  is  considered  to  be 
continuous,  the  constitutive  relationships  arc  available  to  describe  t^.  In  the 
present  formulation,  each  phase  is  assumed  to  be  continuous  and  Newtonian. 
In  ciimpact  tensorial  form,  the  average  stress  vector  for  each  phase  is 


r  ^(HkUkh 

Hk  Tk  =  Pkl - X - 

OXi 


+ 


^{SkiJk)j  ^  2  c  ^{SkUkh 


(/,  j,  n  =  1,2) 


In  x-dircction,  the  average  stress  vector  can  be  expressed  as 


(13) 
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dif'kUk)  2  , 

fik  T^kx  =  l^ik - ^ -  -  T  ^k{  -  + 

(lx  3  (7a' 


^gkVk 

dy 


+  ^^^k  -^^]J  + 

dy 

^  ^kxn  +  '^i^^kf^k^k)  (14) 


If  k^k  i''  assumed  to  be  constant,  it  can  be  shown  that 

^  ■  ^kxo  =  T  ‘  ^  ^ 

Since  the  solid  phase  is  assumed  nondeformable,  it  is  free  of  internal  shear  stress 
(t,  ^  0,  -  0)  and  translates  at  a  prescribed  veloeity,  l<y.  With  this  as¬ 

sumption  and  substituting  Eqs.  (12),  (14),  and  (15)  into  Eq.  (1 1),  it  follows  that 


<^(f>  u) 
Pt 


f  V  -  (p  V  t/)  =  V  -  [///  V(g/  W/)] 


-  V  -  If),  {Vs  -  y)  {us  -  u)  +  Pi  {Vf  -  V)  (m/  -  m)]  -  —  y  p  +  Fj, 

ox 


{  16) 


where 


4 


> 


P  =  SsPs  +  Si  Pi  -  y  /^/^■(;?/  ^i) 


( 17) 


From  Eq.  (3),  it  can  be  shown  that 


Siu,  =  p7  “  -  -p;  Ss^s 


(18) 
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If  it  is  assumed  that  the  phase  densities  are  eonstants  and  the  viscous  stresses 
resulting  from  local  density  gradients  arc  negligible  (V(p/p/)  =  0),  then  from 
Eq.  (IX),  Eq.  (16)  becomes 


//) 

dt 


+  V.(p  1'//)  =  V.(p/  -£-v,/) 


-  V.[p,(E,  -  V){u,  -  u)  +  p,(V,  -  \^{u,  -  JV)]  - 


p  Bx  Fx 


(  19) 


Invoking  the  following  identities 


-  K=fiK 


and 


y  =  fs 


it  can  be  shown  that 


-  u)  +  Pi  (V,  -  y)(u,  -  u)  =  pfji  (20) 


where 


- 


represents  the  relative  phase  velocity.  Substituting  Eq.  (20)  into  Eq.  (19),  it 
follows  that 
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<Xf)  u) 


dt 


+  V  •{p  Vu)  = 


V.(w-^V«) 


pfsfl  KUr  - 


dp 


dx 


p  Bx  +  Fx 


(21  ) 


To  define  the  phase  interaction  force,  Fx,  it  is  necessary  to  consider  the 
multiphase  region  morphology.  For  a  wide  region  of  multi-eonstitutc  solid-li¬ 
quid  phase  change  system,  the  multiphase  region  is  characterized  by  a  fine  per¬ 
meable  solid  matrix.  The  solid  matrix  is  stationary  or  undergoes  free  body 
translation.  Thus  the  liquid  phase  flow  through  the  mushy  region  is  analogy  to 
flow  through  a  porous  media.  Therefore,  implementation  of  Darcy's  law  to 
prescribe  the  phase  interaction  force,  Fx,  is  appropriate.  Thus 


Pi 


imr) 


(22) 


where  Kx  represents  the  component  of  anisotropic  permeability  which  influences 
x-dircction  momentum  transport.  Since 


it  can  be  shown  that 


(23) 

The  second  term  on  the  right-hand  side  of  Eq.  (21)  represents  inertial  forces  in¬ 
duced  as  a  result  of  variations  in  phase  velocities.  This  inertial  force  only  ap¬ 
pears  in  the  multiphase  region,  where  permeabilities  are  extremely  small  and  the 
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inertial  contribution  is  negligible  compared  to  the  Darcian  damping  force. 
Hence  the  x-momentum  retiuccs  to 


<'*(/>  u) 


f  V  -  (p  I'  u) 


V  ■  (/i/  V  w) 


dx 


+  P 


(24) 


This  final  form  of  the  equation  represents  conservation  of  momentum  in  the 
\-direction  for  the  .solid-liquid  phase  change  system.  In  the  solid  region,  the 
permeability  is  zero,  thus  the  velocity  is  the  solid  phase  velocity.  In  the  liquid 
region,  the  pcrmcabilit}  ’s  infinite,  thus  the  Darcian  damping  force  will  disap¬ 
pear.  The  form  o*  :ie  equation  is  the  same  as  the  normal  momentum  equation 
fi)r  a  single  phase  fluid. 


C.  CONSERVATION  OF  ENERGY 


Conservation  of  energy  for  a  two-dimensional,  solid-liquid  phase  change 
system  can  be  expressed  by  the  following  equation: 

(P,  K  +  P,  A,)  +  V  .  (p,  V,  h,  +  p,  y,  AJ  =  V  .  (AV7)  (  25  ) 

where  local  thermodynamic  equilibrium  has  been  assumed  (T*  =  T)  and  the 
mixture  conductivity  is  defined  by 

^  =  8iki  +  gsf<s  (26) 
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The  advcctivc  term  may  be  decomposed  into  contributions  owing  to  the  mean 
mixture  motion  and  the  relative  phase  motion.  Thus 


P/JVb  +  PsV,h,  =  pVh  +  p,{Vi-  mhi-h)  +  pAVs'-  yWh~h) 

(27) 

Where  the  mixture  enthalpy  is 

P/ 

It  is  noted  that  the  flux  owing  to  the  relative  phase  motion  only  has  a  contrib¬ 
ution  in  the  mushy  region.  Substituting  Eqs.(27)  and  (2S)  into  £^.(2.^),  it  follows 
that 

+  V.(pf7z)  =  V.(A:V71 

-  V  -  [p,(  Vi  -  m,  -h)  ^  n  -  n(^  -  m  (29) 

Simplifying  the  term 

p/( ("/  -  n(^  -  /t)  +  Ps( K -  ms  -h)  =  pfjy-  v^h,  -  h,)  oo) 

it  follows  that 

^{ph) 

+  VApVh)  =  V.(/cV71  -  ^-ipfsiV-  mi-h,)-]  (31) 

In  the  present  formulation,  the  enthalpy  of  phase  k  is  defined  as 
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(32) 

•'0 

whore  Ck  represents  an  effective  specific  heat  of  phase  k.  Substituting  the  iden¬ 
tity 


(33) 


into  l:q.(3l),  it  can  be  shown  that 

+  V.(pf7;)  =  V.(-f  V/i)  + V.[-f  V(/i,-/0] 

(7  ‘'i  *'5 


(34) 


This  is  the  final  form  of  the  equation  which  represents  the  conservation  of  en¬ 
ergy  in  the  x-momentum  for  the  solid-liquid  phase  change  system. 
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